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Abstract 

We consider the three-dimensional it J model defined on a simple cubic lattice and study its 
behavior close to the multicritical Nishimori point where the paramagnetic-ferromagnetic, the 
paramagnetic-glassy, and the ferromagnetic-glassy transition lines meet in the T-p phase diagram 
(p characterizes the disorder distribution and gives the fraction of ferromagnetic bonds). For 
this purpose we perform Monte Carlo simulations on cubic lattices of size L < 32 and a finite- 
size scaling analysis of the numerical results. The magnetic-glassy multicritical point is found 
at p* = 0.76820(4), along the Nishimori line given by 2p — 1 = Tanh(J/T). We determine the 
renormalization-group dimensions of the operators that control the renormalization-group flow 
close to the multicritical point, yi = 1.02(5), y2 = 0.61(2), and the susceptibility exponent r] = 
—0.114(3). The temperature and crossover exponents are = l/y2 = 1.64(5) and (f) = yi/y2 = 
1.67(10), respectively. We also investigate the model-A dynamics, obtaining the dynamic critical 
exponent z = 5.0(5). 

PACS numbers: 75.10.Nr, 64.60.Kw, 75.40.-s, 05.10.Ln 
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FIG. 1: Phase diagram of the three-dimensional itJ Ising model in the T-p plane. The phase 
diagram is symmetric for p ^ 1 — P- 

I. INTRODUCTION 

The ±J Ising model provides an interesting theoretical laboratory to study the effects 
of quenched random disorder and frustration in Ising systems. It is defined by the lattice 
Hamiltonian 

^ = - Jxy(^x(^y, (1) 

where a^; = ±1, the sum is over the nearest-neighbor sites of a simple cubic lattice, and the 
exchange interactions Jxy are uncorrelated quenched random variables, taking values ±J 
with probability distribution 

PiJxy) = pSiJxy - J) + (1 - p)6iJxy + J). (2) 

In the following we set J = 1 without loss of generality. For p = 1 we recover the standard 
ferromagnetic Ising model, while for p = 1/2 we obtain the bimodal Ising spin-glass model. 
The ±J Ising model is a simplified model^ for disordered spin systems showing glassy be- 
havior in some region of their phase diagram, such as Fei-^^MnajTiOs and Eui_a;Baa;Mn03, 
see, e.g., Refs. 2,3,4. The random nature of the short-ranged interactions is mimicked by 
nearest-neighbor random bonds. 
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The T-p phase diagram of the three-dimensional ± J Ising model is sketched in Fig. 1 for 
1 > p > 1/2 (it is symmetric for p 1— p)- The high-temperature phase is paramagnetic for 
any p. The nature of the low-temperature phase depends on the value of p: it is ferromagnetic 
for small values of 1 — p, while it is glassy with vanishing magnetization for sufficiently large 
values of 1 — p. The paramagnetic and low-temperature ferromagnetic and glassy phases are 
separated by different transition lines, which meet at a magnetic-glassy multicritical point 
(MGP) located at p*,T* and usually called Nishimori point. 

The paramagnetic-ferromagnetic (PF) transition line starts from the Ising transition at 
p — 1 and extends up to the MGP a,t p — p*. For p — 1 the transition belongs to the Ising 
universality class, while for any 1 > p > p* it belongs to the randomly-dilute Ising (RDIs) 
universality class, ^'^ characterized by the magnetic critical exponents^'^ Vf = 0.683(2) and 
rjf = 0.036(1). The Ising transition at p = 1 is a multicritical point and, close to it, for 
< 1 — p <^ 1, one observes multicritical behavior^'^'^° with crossover exponent = ais, 
where^^ ais — 0.1096(5) is the Ising specific-heat exponent. The paramagnetic-glassy (PG) 
transition hne starts from the MGP and extends up to p = 1/2. A reasonable hypothesis is 
that the critical behavior is independent of p along the PG line, i.e. that a nonzero average 
value [Jxy] of the bond variables is irrelevant at the glass transition, as found in mean-field 
models. Assuming this scenario, for any 1 — p* < p < p* the PG transition belongs to the 
same universality class as that of the bimodal Ising spin glass model at p = 1/2. Its critical 
behavior has been widely investigated (see, e.g., Refs. 13,14 and references therein) and it 
is characterized by the overlap exponents Ug 2.4 and rjg ^ —0.4. 

As argued in Refs. 15,16,17, the MGP is located along the so-called Nishimori line^^'^^ 
(N-line) defined by the relation 

V = Tanh/3 = 2p - 1, (3) 
where /? = 1/T, which allows us to define a Nishimori temperature 

In[P) 2 1-p 

for each value of p. The ± J Ising model along the N-line presents several interesting prop- 
erties. The internal energy has been computed exactly along the N-line:^® 

^iv(p) = :^[(^)T.(p)] = 6p-3, (5) 
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where the angular parentheses and the brackets refer respectively to the thermal average 
and to the quenched average over the bond couplings {Jxy}- Along the N-line several other 
remarkable relations hold, such as^^ 

Ux)] = [{Axn (6) 
where Ax is an arbitrary product of spin variables ax, and also^^ 

G2i+lix) = G2^+2{^) i^l,2,.... (7) 

where Gk{x) = [{o-qo-x)'']. As a consequence of Eq. (7), the magnetic correlation function 
Gi{x) and the overlap correlation function G2{x) are equal along the N-line. The N-line 
separates the regions where magnetic and glassy fluctuations dominate. Arguments based 
on local gauge invariance^^'^^'^^ show that the MGP must be located along the N-line, so that 
T* — Tn{p*). At the MGP, magnetic and glassy fluctuations become critical simultaneously. 
At fixed p an important inequality holds:^^'^^ 

|[(^7,a,-)T]|<[|(a,a,K(,)|], (8) 

where the subscripts indicate the temperature of the thermal average. This relation shows 
that ferromagnetism can exist only in the region p > p* and that the system is maximally 
magnetized along the N-line. Ref. 20 (see also Refs. 19,21) also reports an argument that 
indicates that the ferromagnetic-glassy (FG) transition fine coincides with the fine p — p*, 
from T = T* to T = 0. This conjecture is contradicted by recent results for the two- 
dimensional ±J modeP^'^^'^^ and for three-dimensional random-plaquette gauge model,^^ 
which is the dual of the ±J model. Violations are in any case quite small. We mention 
that a mixed low-temperature phase, in which ferromagnetism and glass order coexist, is 
found in mean-field models^^ such as the infinite-range Sherrington-Kirkpatrick model.^^ Its 
presence has been confirmed in the ±J Ising model defined on Bethe lattices. However, 
there is no evidence of this mixed phase in the ± J Ising model on a cubic lattice^^ and in 
related models. Nevertheless, the existence of such a mixed phase is still an open problem, 
as discussed in Ref. 27. 

In this paper we consider the ±J model and perform Monte Carlo (MC) simulations 
along the N-line close to the MGP. By performing a finite-size scaling (FSS) analysis, we 
locate the multicritical point along the N-line, finding p* — 0.76820(4). We determine the 
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renormalization-group (RG) dimensions yi and 1/2 of the relevant operators that control the 
RG flow close to the MGP and the exponent rj that gives the critical behavior of the magnetic 
and of the overlap susceptibility. We obtain yi = 1.02(5), y2 — 0.61(2), and 77 = —0.114(3). 
The temperature and crossover exponents are u — l/y2 — 1.64(5) and (f) — y2/yi — 1.67(10) 
respectively. We also use our numerical results to estimate the dynamic critical exponent z 
that characterizes the model-A dynamics^° at the MGP, i.e. a relaxational dynamics without 
conserved order parameters. We obtain z = 5.0(5). Our results significantly improve those 
obtained in previous -y^^orks.^^'^^'^^'^^'^^'^^ 

The paper is organized as follows. In Sec. II we summarize the theoretical results we 
need in our numerical analysis. In Sec. Ill we report our numerical results. We estimate 
the position of the MGP and the critical exponents yi, ?/2, and 77 in Sec. Ill A, while in 
Sec. IllB we give an estimate of the exponent z for the Metropolis dynamics we use, which 
is a specific example of a relaxational dynamics without order parameters (the so-called 
model-A dynamics). In Sec. IV we summarize our results. In the Appendix we report some 
notations. 

II. SUMMARY OF THEORETICAL RESULTS 

In the absence of external fields, the critical behavior at the MGP is characterized by 
two relevant RG operators. The singular part of the free energy averaged over disorder in a 
volume of size L can be written as 

F,,^,{T,p,L) = L'"'f{u,Ly\u2Ly\{u,Ly^}), t>3, (9) 

where yi > 1/2 > 0, yi < for i > 3, Ui are the corresponding scaling fields, and Ui = U2 = 
at the MGP. In the infinite-volume limit and neglecting subleading corrections, we have 

F,i„g(T,p) = \u2\''^y'f±{u,\u2r'^), <P = 1/1/1/2 > 1, (10) 

where the functions f±{x) apply to the parameter regions in which ±1^2 > 0. Close to the 
MGP, all transition fines correspond to constant values of the product Ui\u2\~'^ and thus, 
since > 1, they are tangent to the line Ui = 0. 

The scaling fields Ui are analytic functions of the model parameters T and p. Using 
symmetry arguments, Refs. 16,17 showed that one scaling axis is along the N-line, i.e. that 
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the N-line is either tangent to the hne Mi = or to U2 = 0. Since the N-hne cannot be 
tangent to the transition hnes at the MGP and these hnes are tangent to ui — 0, the first 
possibihty is excluded. Thus, close to the MGP the N-line corresponds to U2 — 0. Thus, we 
identifyi^'i^ 

U2 = v-2p+l. (11) 

As for the scaling axis Ui — 0, e = 6 — d expansion calculations predict it^^ to be parallel to 
the T axis. The extension of this result to d — 3 suggests 

ui^p-p*. (12) 

Note that, if Eq. (12) holds, only the scaling field U2 depends on the temperature T. We 
may then identify u — l/y2 as the critical exponent associated with the temperature, and 
rewrite Eq. (10) as 

Fsin,(T,p) = |i|'^V±(^in, (13) 

where t = {T — T*)/T*, g = p — p*, and is the crossover exponent. 

These results give rise to the following predictions for the ESS behavior around T*, p*. 
Let us consider a RG invariant quantity R, such as = ^/ L, U4, U22, which are defined in 
the Appendix, and its derivative R' with respect to /3 = l/T. In general, in the ESS limit 
R obeys the scaling law 

R = n(uiLy\u2Ly\{uiLy'}), i > 3. (14) 

Neglecting the scaling corrections, that is terms vanishing in the limit L — > 00, close to the 
MGP we expect 

R = R* + buUiLy^ + h2iU2U'^ + . . . . (15) 

which is valid as long as uiW'^ is small. Along the N-line, the scaling field U2 vanishes, so 
that we can write 

Rn^R* + bnuiiy^ + ..., (16) 

where the subscript N indicates that R is restricted to the N-linc. Let us now consider the 
derivative of R with respect to p. Differentiating Eq. (15), we obtain 

R' = biiu[Ly' + 62i«2^^' + • • • (17) 
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If Eq. (12) holds, then u[ = 0, so that 



R' = 621 «2^^' + • • • 



(18) 



This result gives us a method to verify the conjecture of Ref. 17: once yi has been determined 
from the scaling behavior of a RG invariant ratio close to the MGP, it is enough to check 
the scaling behavior of R'. If R' scales as with x < yi, the conjecture is confirmed and x 
provides an estimate of y2- 

Finally, we consider the magnetic susceptibihty. Along the N-line it behaves as 



Note that there is only one rj exponent which characterizes the critical behavior of both 
the magnetic and overlap correlation functions, since they are equal along the N-line, see 



In the following we present a FSS analysis of high-statistics MC data along the N-hne 
close to the MGP. We performed MC simulations for lattice sizes L = 8, 12, 16, 24, 32, taking 
periodic boundary conditions. Wc used a standard Metropolis algorithm and multispin 
coding (details can be found in Ref. 6) . Most of the simulations correspond to values of p 
in the range 0.7680 < p < 0.7685, i.e. very close to the MGP, which, as we show below, is 
located at p* — 0.76820(4): typically, we considered 6 values of p in this range for each value 
of L. To obtain small statistical errors, we generated a large number of samples: 2 x 10^ for 
L < 16, 10^ for L = 24, and 4 x IC^ for L = 32. Because of the long equilibration times, 
for each sample wc performed a large number of Metropolis sweeps; for L = 16, 24, 32, 
the number of sweeps is 10^, 8 x 10^, and 5 x 10^, respectively. To guarantee equilibration, 
typically 30% of the data were discarded (but, for L — 32, we discarded 50% of the data). All 
MC data are available on request. Below we report the results of the analyses: in Sec. Ill A 
we consider the static exponents, while in Sec. Ill B we focus on the dynamics. 



XN^eL^-^{l + e,u,Ly' + ■■■). 



(19) 



Eq. (7). 



III. RESULTS 
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FIG. 2: MC data of = i/L, U4, and Ud vs p. 
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P* yi u^2 Rl 






0.88 


0.59910(2) 1.02(5) 0.3180(3) 0.5648(3) 




U22, R^, 


Ud 1.43 


0.59902(2) 1.02(4) 0.3189(3) 0.5640(4) [/* 


= 1.2137(3) 


U22, R^, 


C/4 0.62 


0.59914(2) 1.01(5) 0.3178(3) 0.5656(4) 


= 1.5302(6) 



TABLE I: Results of combined fits. The first fit uses all data with L > 8, the last two fits only 
those with L > 12. DOF is the number of degrees of freedom of the fit. 

A. Static exponents 

MC estimates of the RG invariant quantities i?^, C/4, and Ua along the N-line are shown in 
Fig. 2. There is clearly a crossing point at p ^ 0.7682, which provides a first rough estimate 
of the location of the MGP point. In order to estimate precisely p*, T*, and yi we fit the 
renormalized couplings R close to the MGP to 

R = R* + a{f3- f3*)Ly\ (20) 

keeping R*, and yi as free parameters. Note that this functional form relies on the 
property that U2 = along the N-line. Otherwise, an additional term of the form (/? — /3*)L^2 
should be added. We also neglect scaling corrections that behave as cL^-^ with y-^ < 0. Indeed, 
since we only have data in a limited range of values of L, we are not able to include reliably 
a correction of this type. 

Fits that involve R^ and U22 have an acceptable even if we include all data with L > 8: 
there is no evidence of scaling corrections. On the other hand, in fits of C/4 or the data 
with L = 8 must be discarded to obtain a good x^- To obtain more accurate estimates, we 
have performed combined fits in which several RG invariant quantities are fitted together. 
The results are reported in Table 1. The dependence on the observables used in the fit is 
reasonably small and allows us to estimate 

= 0.5991(1), (21) 
yi = 1.02(5). (22) 

The errors take into account the variation of the estimates with the different observables 
used in the fits (note that statistical errors are much smaller). Since scaling corrections are 
expected to differ in the different observables, this should allow us to take indirectly into 
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TABLE II: Estimates of x. We report results obtained by analyzing simultaneously two different 
quantities and including only data satisfying L > Lmin- DOF is the number of degrees of freedom 
of the fit. 

account the scaling corrections. We have then T* — 1/P* — 1.6692(3), and, by using Eq. (3), 

p* = 0.76820(4). (23) 

In Table I we also report estimates of the critical value of the RG renormalized couplings. 
Note that U22 ~ 0.318, which is significantly higher than the corresponding result for 
the RDIs universahty class, U22 = 0.1479(6).^ This indicates^'' that the violations of self- 
averaging are much stronger at the MGP than along the PF transition line, as of course 
should be expected. 

We consider now the derivative R' of the RG invariant quantities with respect to p. They 
have been determined by considering the connected correlations of R and of the Hamiltonian. 

At the critical point, R' is expected to behave as for large L, where x = y2, if the argument 
of Ref. 17 holds; otherwise, one should have x = yi. In order to determine x, we fit InR' to 

InR' ^a + x\nL + b{(3 - (3*)Ly\ (24) 

keeping yi fixed to yi — 1.02(5). To avoid fixing /3c we perform combined fits in which 
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one derivative R'^ and one RG coupling R2 are fitted together. The resuhs are reported in 
Table II. The of the fit is always good except when we use Lmin = 8 and C/4. If we do 
not consider the corresponding results, all estimates of x are close to 0.61. Analyses of R'^ 
are apparently stable with Lmin, while those of C/4 show a slight upward trend. A reasonable 
final estimate \s x — 0.61(2), which takes into account all results with their error bars. This 
result is significantly different from yi and thus confirms the argument of Ref. 17. Since 
X < Hi, X should be identified with 1/2. Therefore, we obtain the estimates 

7/2 = 0.61(2), i/= — = 1.64(5). (25) 

The crossover exponent is therefore 

0=^ = 1.67(10). (26) 

The same analysis used to estimate 1/2 can be employed to determine 77. Instead of x, we 
consider the ratio Z = x/C^, which has smaller statistical errors. Since Z ~ for L — > 00 
at the critical point, we fit the MC data to 

In Z = a - 77 In L + 6(/3 - /3*)L^i . 

As before, we fix yi and perform combined fits of In Z with a RG invariant coupling, consid- 
ering only data satisfying L > Lmin- Fits of Z and R^ give 77 = —0.1155(6) and —0.1154(9) 
for Lmin = 8, 12; if we use U22 instead of R^, we obtain r] — —0.1134(7) and —0.1131(9) 
for Lmin = 8, 12. The Lmin dependence is small and results change only slightly with the 
observable. We take as our final estimate 

77 = -0.114(3). (27) 

Our FSS results significantly improve earlier results. Ref. 33 reports the computation and 
analysis of the 34th-order high-temperature (HT) series of some susceptibilities 

Xm,n=l^J2^{SiS,rr (28) 
ij 

along the N-line, obtaining p* = 0.7656(20), yi = 1.18(11), (p = yi/y2 = 1.85(14), t] = 
—0.10(2). These estimates are substantially consistent with ours. As a further check, we 
reanalize the 34th-order HT series reported in Ref. 33, by biasing the value of the critical 
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FIG. 3: Estimates of the exponential autocorrelation time r at the MGP vs L. 

point with the MC estimate (21). Using biased first-order integral approximants, see, e.g., 
Ref. 11 for details, we obtain (2 — 'r])ly\ = 2.08(7) from the series of xii, (1 — '^v)/yi = 
1.25(17) from the series of X22, 3/?/i = 3.03(14) from the series of the ratio X11/X22, and 
{2—r]—y2)/yi = 2.70(9) from f(9x2i/<9f , from which we can derive the estimates yi = 0.99(5), 
= yi/y2 = 1-6(3), and 77 = —0.1(1), which are in good agreement with our FSS results. 

Other results can be found in Refs. 31,32,35; they are apparently less precise and not con- 
sistent with ours within the reported errors. For example, we mention the recent estimates 
p* = 0.7673(3) obtained by off-equilibrium MC simulations^^ and p* ~ 0.622 obtained by a 
RG study. Note that estimate (23) and the conjecture^^ of Refs. 39,40 allow us to find the 
location of the multicritical point that occurs in the three-dimensional random-plaquette 
gauge model. We obtain Pg^uge = 0.9650(1), which is in agreement with, though much more 
precise than, the result of Ref. 41, p* „e = 0.967(4). 



B. Model- A dynamic exponent z 

Finally, we present some results on the dynamic behavior of the Metropolis algorithm, 
which represents a particular implementation of a relaxational dynamics without conserved 
order parameters (model-A dynamics). Note that at the MGP there is only one dynamic 
exponent z characterizing the relaxation of both the magnetic and the glassy critical modes, 
since their autocorrelation functions are strictly equal along the N-line.^^ In Fig. 3 we show 
estimates of the exponential autocorrelation time r at the MGP as extracted from the 
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connected autocorrelation function of the magnetic susceptibility 

G^{t,-h)^[{x{h)x{t2))c]. (29) 

For large L and T = T*, r is expected to scale as L^, where z is the dynamic critical 

exponent. A hnear fit of the MC results to Inr = a + 61nL gives the estimate z = 5.0(5), 

which is significantly larger than the value at the PF transition line z = 2.35(2).'^^ Instead 

this estimate is close to the value of z obtained for the bimodal Ising spin-glass model, 
^ = 5.7(2).43 

We also determine the exponent A which describes the nonequilibrium relaxation of the 
magnetization at Tc from a starting configuration in which all spins are parallel. Asymp- 
totically, for f — > oo, one expects 

M(t) ~ ^-^ A = (30) 

iz 

see, e.g., Ref. 36 and references therein. Our results lead to the estimate A = 0.09(1), 
which is perfectly consistent with the estimate^^ A = 0.090(3) obtained in off-equihbrium 
MC simulations. 

IV. CONCLUSIONS 

In this paper we have considered the critical behavior close to the MGP which is present 
in the phase diagram of the ± J model. Our main results are the following: 

(i) We have obtained an accurate estimate of the location of the MGP: -p* = 0.76820(4), 
/3* = 0.5991(1). It is worth observing that our estimate of p* is very close to the 
result^^ Pc = 0.778(5) for the location of the FG transition at T = and satisfies the 
rigorous inequality Pc > V* which follows from Eq. (8). Our results show therefore 
that, even if the conjecture^^'^^ that the FG transition hue does not depend on the 
temperature is not true, deviations are quite small. 

(ii) We have verified the conjecture of Ref. 17: the scaling field U\ associated with the RG 
operator with the largest RG dimension does not depend on the temperature. 

(iii) We have determined the critical exponents y\ — 1.02(5), y2 — 0.61(2), v = l/y2 — 
1.64(5), = yi/y2 = 1.67(10), and 77 = -0.114(3). 
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(iv) We have determined the dynamic critical exponent z associated with the model-A 
dynamics, obtaining z — 5.0(5). 

Our resuhs are significantly more precise than those obtained in previous works. ^^'^^'^^'^^'^^'^^ 
They can be used to explain experiments on materials containing both ferromagnetic and 
antiferromagnetic ions. An example is Fej;Mni_3;Ti03, which shows Ising behavior for x = 1 
and X = 0, a PG transition for 0.38 < x < 0.58 and a PF transition for < x < 0.38 and 
0.58 <x< l.^'''^^ The MGP should be located at x fa 0.38 and at x fa 0.58. Close to these 
values our results apply. 

Acknowledgments 

The MC simulations have been done at the Computer Laboratory of the Physics Depart- 
ment of Pisa. 

APPENDIX A: NOTATIONS 



where the angular parentheses and the brackets indicate respectively the thermal average 
and the quenched average over J^y^ the magnetic and overlap correlation functions are given 
respectively by Gxix) and (a;). Along the N-line, cf. Eq. (3), G\{x) = (x). 
We define the magnetic susceptibility x = Sa;G^i(^) the correlation length ^ 



where gmin = (27r/L, 0,0), q = 2 sing/2, and Gi{q) is the Fourier transform of Gi{x). We 
also consider quantities that are invariant under RG transformations in the critical limit. 
Beside the ratio 



Setting 



(Al) 




(A2) 



(A3) 



we consider the quartic cumulants t/4, t/22, and Ud defined by 



t/4 = 



(A4) 



2' 
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^22 = 



where 

(A5) 

X 

Analogous quantities C/4, U22: and can be defined by using the overlap variable 
= al^^a^\ where the superscripts indicate two independent configurations for given 

disorder. Using Eq. (6), one can easily check that along the N-line = and C/4 = C/4. 

This implies that also their fixed-point values are the same at the MGP. 
Finally, we consider the derivatives 



which can be computed by measuring appropriate expectation values at fixed (5 and p. 
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